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Estimation of process variables 



DESCRIPTION 

Field of the Invention 

The invention relates to the field of process control. It relates to a method, 
computer program and data processing system for the estimation of process 
variables as described in the preamble of claims 1 , 8 and 9, respectively. 

Background of the Invention 

Mathematical modelling is an essential tool in modern industry. It makes 
possible a whole range of procedures for process control, planning, 
scheduling, optimisation, monitoring etc., where processes considered are 
technical or physical processes such as industrial production processes, power 
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plants, mechatronic devices etc. The success of the above procedures hinges 
on the ability of the mathematical models to represent the physical reality at a 
needed degree of accuracy. The mathematical models incorporate a number of 
parameters. Some of these parameters have a physical meaning, e.g. mass, 
5 flow, temperature, force. Other parameters are coefficients of given 
parameterised functional dependencies among physical properties of the 
process. For example, in a turbomachine, i.e. a compressor or a turbine, such 
dependencies are expressed by polynomials describing relationships such as 

io mass flow =fl (vane position, pressure ratio, rotor speed, etc) 

efficiency=f2(temperature, mass flows, pressures, etc) 

where fl and fl are static non-linear functions of the bracketed variables, 
which represent measured or estimated process states in the machine. 

is Graphical representations of these functions are so called "compressor maps". 
Coefficients of polynomials replacing these functions are usually computed via 
least square fitting of experimental data. Usually only manufacturers have 
enough data to carry out the fitting at all relevant conditions. Moreover, the 
process of fitting is costly and cumbersome. During operation, not in a test 

20 set-up but at an arbitrary installation, it is difficult to measure accurately all 
the values entering the relationships fl and f2 and therefore to deduce a 
suitable approximation. 

In order to obtain estimates of the relationships on-line, i.e. during process 
25 operation, estimation techniques have been. developed: The "Kalman Filter" (KF) 
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for linear systems emerged in the 1960's and has been developed into a 
mature technology for reconstruction of the state of a given linear plant. The 
"Extended Kalman Filter" (EKF) is used when the KF formulae are applied to 
estimation of the internal states of nonlinear plants. In this case, linearisations 
of the true system equations, at current estimates, are used. If parameters 
entering the plant are unknown then the "State Augmented Extended Kalman 
Filter" (SAEKF) is used. The SAEKF provides estimates of the polynomial 
coefficients representing unknown relationships within the mathematical 
process model, such as the functions fl and f2 mentioned above. For example, 
f2 could be represented by a third-order polynomial 
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where ko ... k 3 are vectors with polynomial coefficients to be estimated, and xi 
... Xn are process states representing vane position, pressure, speed, etc. The 
complete estimation of a physical property of the process such as the 
efficiency e is done according to the block diagram of figure 1 . The polynomial 
relation shown above is represented by the right hand side of the block 
diagram, comprising multiplication and addition blocks. A state vector 
according to the SAEKF comprises the coefficient vectors k 0 ... k 3 as extended 
state variables. In the block diagram, these states are represented by outputs 
of integrators having a constant input of zero. The SAEKF algorithm modifies 
the states, providing updated values for the coefficients ko ... k3, which are 
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then combined with the states xi ... x„to compute an estimate of the efficiency 
e. 

However, because the coefficients ko ... k 3 lack physical meaning, their values 
5 are very sensitive to variations in the data and are not adequate for 
extrapolation to conditions other than the ones they were generated from. 
Furthermore, the polynomial parameterisation restrains the degrees of 
freedom, which might be crucial for an adequate description of the physical 
phenomena under investigation. For instance, the real physical relationships 
10 might incorporate discontinuities, which would not be adequately represented 
if the experimental data is represented by polynomial fitting. 

Description of the Invention 

15 

It is therefore an object of the invention to create a method, computer program 
and data processing system for the estimation of process variables of the type 
mentioned initially, which overcomes the disadvantages mentioned above. 

20 These objects are achieved by a method, computer program and data 
processing system for the estimation of process variables according to the 
claims 1 , 8 and 9. 

In the inventive method for estimating a value of a vector of variables p in a 
25 mathematical model representing a physical process, the vector p represents 
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one or more properties of the process where these properties are 
representable by a function of the state vector, and are estimated as an 
augmented state by the SAEKF. 

5 That is, the process properties such as mass flow or efficiency themselves are 
estimated, and not polynomial coefficients for computing the variables from 
the state! The surprising result is that this approach works and that the 
abovementioned disadvantages are overcome. 

10 In a preferred variant of the invention, the physical process comprises a 
turbomachine, and the vector of variables p represents at least one of an 
efficiency or a mass flow rate of the turbomachine. 

In a further preferred variant of the invention, the physical process comprises a 
15 heat exchanger, and the vector of variables p represents at least one heat 
transfer coefficient of the heat exchanger. 

In a further preferred variant of the invention, the physical process comprises a 
mating gear transmission, and the vector of variables p represents a backlash 
20 and spring function. 

The computer program for the estimation of process variables according to the 
invention is loadable into an internal memory of a digital computer, and 
comprises computer program code means to make, when said computer 
25 program code means is loaded in the computer, the computer execute the 
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method according to the invention. In a preferred embodiment of the 
invention, a computer program product comprises a computer readable 
medium, having the computer program code means recorded thereon. 

A data processing system for the estimation of process variables comprises 
means for carrying out the steps of the method according to the invention. In a 
preferred embodiment of the invention, the data processing system is an 
apparatus comprising a data processor, a memory coupled to the processor 
and computer program code means stored in said memory, where said 
computer program code means, when executed by the processor, causes the 
method according to the invention to be executed. 

Further preferred embodiments are evident from the dependent patent claims. 

Brief Description of the Drawings 

The subject matter of the invention will be explained in more detail in the 
following text with reference to preferred exemplary embodiments which are 
illustrated in the attached drawings, in which: 

Figure 1 shows a block diagram for the computation of a process property 

according to the state of the art; 
Figure 2 shows a block diagram for the computation of a process property 

according to the invention; 
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Figure 3 shows a block diagram representing the computation of process 

variables in a flow process; 
Figure 4 schematically shows process variables of a heat transfer process; 
Figure 5 schematically shows a mechanical system with backlash; 
Figure 6 shows a backlash torque tk with memory effect; 
Figure 7 shows a polynomial approximation to a backlash and spring 

function; and 

Figure 8 shows a continuous model approximation of friction torque. 

The reference symbols used in the drawings, and their meanings, are listed in 
summary form in the list of reference symbols. In principle, identical parts are 
provided with the same reference symbols in the figures. 



Detailed Description of Preferred Embodiments 

The invention shall first be explained in general terms as a special embodiment 
of a State Augmented Extended Kalman Filter (SAEKF), and then as applied to 
specific examples, i.e. a turbomachine, heat exchanger and gear backlash. 

The well-known Kalman Filter estimates a system state of a dynamic system 
that is representable by state-space representation 
x = f(x,u) 

with state vector x and input vector u. In the SAEKF, the state is augmented by 
variables p to be estimated, and the state-space model underlying the SAEKF is 
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where x = f(x,u,p) represents a known dependency of the change x in system 
state from the state vector x, the measured values u and the vector of variables 
p, and ^represents a vector of noise disturbances. The original state vector x 
5 shall be referred to as state vector, whereas the vector of variables p is 
denoted as augmented state. The combined vector 
x 

_P„ 

shall be called complete state. Estimates of the states are denoted by xand p. 
In a preferred embodiment of the invention, the variables p to be estimated 
10 correspond to physical properties of the process, i.e. to properties whose value 
have an interpretation in terms of the process, as opposed to e.g. coefficients 
of the polynomial approximations. In the simplest case, the vector p is a scalar. 

The computation of the state estimate along with associated covariance 
15 matrices is done according to a suitable implementation of the known SAEKF 
approach, as shown e.g. in Robert Stengel, "Optimal control and estimation", 
Dover Publications, 1994.; and in C. Bohn, "Recursive Parameter Estimation for 
Nonlinear Continuous Time Systems through Sensitivity Model Based Adaptive 
Filters", PhD Dissertation, University of Bochum, Germany, 2000. The Extended 
20 Kalman Filter is closely related to a wide family of system identification 
methods called Recursive-Prediction-Error-Methods (RPEM). The main 
difference between the methods in this class lies in the computation of the 
gradients of the optimisation criterion, see for example L. Ljung, ""System 
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Identification: Theory for the User", 1999, Chapter 11. Hence, the concept of 
replacing polynomial approximations with one single augmented state also 
applies in the more general setting of RPEM. 

In a preferred embodiment of the invention, the computations according to the 
SAEKF are executed as summarised by the following steps. Numerous other 
implementations of the SAEKF are possible, with different numerical accuracy, 
effort and stability. 

1. State estimate propagation, solved e.g. by numeric integration of the 
complete state equations for the time fin the range t k <f <r* +1 : 



2. Covariance propagation, solved either by numeric integration or by explicit 
solution of the differential equations: 



x = f(x , u, p(t k )) 




xp 



dp 




pp 



3. Predicted output computation: 



K'*+l ) = h {x(.t* + \ )> «('**! )>pQt )) 
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4. Prediction error computation: 

5. Approximate prediction of error covariance matrix computation: 

dh , ,dh T dh n f+ . .dh T 

+ y" p * ( ' w) a* + ^ T ^ ** ,; dp 

6. Computation of Filter gain matrices, where K is the Kalman Filter Gain for 
the state x and L is the Kalman Filter Gain for the augmented state p: 

7. State update computation: 

x(/; +1 ) = x(f; +1 )+Ai(^ +1 ) 

8. Parameter Update computation: 

Apa* t i)= i ('wM'w) 
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9. Covariance update computation: 

^,(0 = P xp (t- M )-K(t t+l )A(t k+l )L r (t k+l ) 

p„ (C ) = p pp (t™ ) - Uf** ) A ih^ (f M ) 

The above steps are performed iteratively, for increasing time t and with 
measured variables u. What makes the SAEKF work, is the connection of the 
variables p and the system states x through the covariances of the noise 
processes. 



Turbomachine application: 

A compressor map typically describes a dependency of the mass flow m of a 
working fluid through a compressor or turbine, and/or an efficiency e of the 
compressor or turbine process. Usually, a pressure ratio, a rotational speed 
and optionally a guide vane angle are measured or estimated, and from these 
the mass flow and efficiency are determined. This is done either manually from 
a graphic representation of the map, or in a control system according to look- 
up tables and computations, or from polynomial approximations to the 
respective functions. The dependency of the efficiency and mass flow on the 
known variables is therefore static and nonlinear. Due to degradation of the 
machine, e.g. by blade erosion, the dependency changes over time, and a 
continuously updated estimate of the dependency is desired. 
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For this estimation, the traditional approach is to model the dependency of 
mass flow on the system state as 

m^m ref C x {w)C 2 {v) 

Cj(w) = a Q + a Y w + a 2 w 2 + a 3 w 3 + • • • 

c 2 ( v ) = fc 0 + v+ V 2+ V 3 + ••• 

where m is a normalised mass flow, w represents rotational speed and v 
represents a guide vane angle. From m and the pressure ratio, the estimate of 
the actual mass flow m is computed. Alternatively, the dependency on the 
pressure ratio dP can be expressed as a third polynomial as illustrated below. 

m = m ref C l (w)C 2 (v)C 3 ( < dP) 

Q (w) = a 0 + a x w + a 2 w 2 + a 3 w 3 + • • • 

c 2 (v) = ^ 0 +v+V 2+ V 3+ -*- 

C 3 (dP) = c 0 + c x dP + c 2 dP 2 + c 3 dP 3 + — 

The SAEKF is applied to estimate the polynomial coefficients ao ... a n and bo ... 
b n . From a current estimate of said coefficients and from current values of w, v 
and pressure ratio, the mass flow m is determined. An analogous procedure is 
applied for the efficiency e. 

According to the invention, it is not necessary to have an explicit 
representation of the dependency of m and e on the system state. It is 
sufficient to know that there exists such a dependency. The dependency can 
be either linear or nonlinear, static or with memory. An example of a nonlinear 
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dependency with memory is a backlash model of a gear transmission as 
represented by Figure 6. 

According to the invention, the variables m and e are considered as 
augmented states of the corresponding compressor or turbine model. The 
SAEKF approach, or a suitable modification, is then used to estimate the 
augmented states. The SAEKF then tracks the mass flow m and efficiency e, 
although they change both relatively fast with the operating state of the 
machine, and relatively slow with a degradation of the machine. 

In a preferred embodiment of the invention, the variables estimated in this way 
are influenced by conditions that change slowly with respect to the main 
process dynamics, e.g. turbine degradation due to erosion and fouling. The 
efficiency estimate under the same operating conditions shall deteriorate over 
time, and its value can therefore be used as a basis for preventive 
maintenance, equipment health monitoring, fault diagnosis, etc. 

The operating state or condition undergoes changes with time constants in the 
range of less than a second to several minutes, whereas the degradation has a 
noticeable effect over weeks and months. The traditional approach would only 
compensate for the degradation, whereas the state changes would be taken 
into account through the computation of m and e from the state through the 
polynomially approximated functions. 
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In another preferred embodiment of the invention, an explicit representation of 
the dependency of m and e is determined from a parametric fitting of the 
estimates provided by the SAEKF and the corresponding measured or 
estimated values of the arguments such as speed, pressure ratio, etc. The 
5 explicit model then can be used in offline "what-if" simulations. 

The representation according to equation (1), as used in the SAEKF is obtained 
in the following way: Figure 3 shows a block diagram representing the 
computation of process variables in a flow process. Arrows entering a block 

10 denote input variables to the block, from which the block computes output 
variables denoted by outgoing arrows. N1 is a flow element which computes a 
flow Fl through e.g. a turbine or compressor and an exit temperature Tl i from 
inlet pressure P0, outlet pressure PI, inlet temperature TO, rotational speed w 
and guide vane angle v. The computation corresponds to nonlinear and static, 

15 i.e. memoryless functions. Analogous relations hold for flow element N2. A 
flow element may also represent a valve, with an input corresponding to valve 
stroke instead of w and v, or simply part of a pipe. 

VI is a storage element, in which a pressure PI and temperature Tlo in the 
20 element are computed from inlet temperature Tl i and from an inflow Fl and 
outflow F2. These computations are based on mass and energy balances in the 
storage element, and represent the dynamic part of the process. In all 
relations, enthalpies may be used instead of or in addition to temperatures. 
Flows are mass flows m or volume flows. 



25 
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Flow and storage elements are daisy-chained as in Figure 3 to represent a 
complete compressor or turbo-group together with preceding and subsequent 
piping. In order to obtain the state-space representation of equation (1), the 
states of the storage elements are collected in the state vector x, and the flows 
5 and optionally the efficiencies are collected in the augmented state vector p. In 
the vector equation x = f(x,u,p), the vector function f represents the modelled 

dynamics of all the storage elements V1,V2,... Such dynamic models for 
turbines and thermohydraulic processes are commonly known. According to 
the invention, the static relations embodied by the flow elements Nl ,N2,... and 
10 leading to the vector of variables p are not modelled explicitly. Instead, the 
variables in p, e.g. the flows, are determined by the SAEKF. The derivatives of 
the variables in p are set to zero, as expressed by equation (1), and their 
values are updated by the SAEKF through the effect of the covariance matrices 

Pxx, Pxp, Ppp- 

15 

Another example of a polynomial approximation is found in Knud Rasmussen, 
"Calculation methods for the physical properties of air used in the calibration 
of microphones", Department of acoustic technology, Technical university of 
Denmark report PL-lib, 1997 (http://www.dat.dtu.dk/docs/PLllb-RAP.PDF) 
20 where the specific heat ratio of humid air is approximated as follows: 

k = ao + ai T + a 2 T 2 + {a 3 +a 4 T + asT 2 )x w + (a G + a 7 T + asT 2 )P + a^i + a 10 P 2 

P SV (T) = earp(6 0 r 2 + 6 1 T+6 2 + 6 3 T- 1 ) 
f(P,T) = c 0 + Cl P + c 2 T 2 
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where ao, ai, a2, a3, a4, as, a6, a7, as, ag, aio, bo, bi, b2, b3, co, ci, C2 are 
polynomial coefficients to be estimated. In the formulas above, k is the specific 
heat ratio, T is temperature, P is pressure, x w is a mole fraction of water vapour 
in air, Psv is the saturation water pressure and H is the relative humidity in 
5 percentage. According to the invention, the polynomial for the specific heat 
capacity k is replaced with one single augmented state and included in a larger 
model, for example in a model of a compressor, a turbine or a heat exchanger. 

Heat exchanger application: 

10 

Figure 4 illustrates a model for a heat transfer process in a heat exchanger 
pipe or pipe section: Hot exhaust gas with an inlet temperature T gi and mass 
flow rate f g transfers energy to a metal pipe wall with temperature T m . The 
energy is further transferred from the pipe to steam with an inlet temperature 
15 Tsi and flow rate f s . The gas temperature is changed to T g and the steam 
temperature to T s . The following equations may be used to model the process, 
where heat transfer coefficients from gas to metal and metal to steam are 
denoted by ctg and a s , respectively: 

t g =<x g V m -T g ) + f g (T gi -T g ) 
20 f m = a g (T g -T m ) + a s (T s -T m ) 
t s =a s (T m -T s ) + f s (T si -T s ) 

The inlet temperatures T si , T gi are assumed to be known, the temperatures T g , 
T m , T s are model states. Alternatively, enthalpy is used as a steam state and 
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steam temperature is computed from steam tables. Some or all of the states 
are measured, with measurement noise vector v. 

The known approach would be to represent otg and a s by polynomials in terms 
5 of the flow rates f g and f s . , i.e. 
a g =a 0 + a x f g +a 2 f g +... 

a, = fi 0 + Af. + /hf? + - 
and to estimate the polynomial coefficients ceo , ai , ... and po , pi with a SAEKF. 
According to the invention, the heat transfer coefficients otg and oc s themselves 
are incorporated as augmented states in the vector of variables p and 
10 estimated by the SAEKF. 

Gear backlash application: 

The differential equations of motion for an exemplary two-mass motor-arm 
5 system in a robot, schematically shown in Figure 5, are as follows. 

&i — #3 
#2 = X4 

Xz = — (u — D(x 3 - X±) — TK—Tf — K r ) 

mi 

X\ = — (D(X 3 - X 4 ) + T K ) 

m.2 

T f = C v x 3 + C c — atan(ax 3 ) 

where xi, x 2 are motor and arm positions, and mi, rri2 are motor and arm 
inertias, respectively, tk is a nonlinear backlash and spring function to be 
identified. K r represents a proportional or P-controller on the motor position, u 
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is an additional torque input used for identification, D is a known damper 
coefficient and x f is a motor friction torque. The model is written on short form 
as follows 

x = f (x, u, p) + v 

where P =[tK]T is the augmented state and the unknown backlash and spring 
function to be estimated. Figure 6 shows a backlash torque xk with memory 
effect as a function of a position difference xl-x2 between the motor M and 
the arm A. x f is the viscous and Coloumb friction torque model. To get the 
Kalman filter to work properly, a continuous approximation of the Coloumb 
friction is needed in order to compute the gradients df/dx at speeds close to 
zero. The friction torque as represented by the above equation for x f is 
illustrated in Figure 8. 

The measurements are either motor velocity only or both motor and arm 
velocity. The two candidates for the measurement vector y equal 

y 2 = [X Z +Wl,X 4 -\-W 2 }' 

where w, w, and w 2 are noise processes. The complete model is illustrated in 
Figure 5. 

One advantage of the State-Augmented filter approach according to the 
invention is the fact that it is not necessary to assume a particular model 
structure for the backlash and spring function to be identified. Instead, tr is 
simply made an augmented state. 



iT 



02/097 



- 19 - 

The augmented states are introduced into the system equations in 
the following way. 

f(x,u,p) 1 
0 J 

5 Traditionally, polynomials have been used to approximate nonlinear functions 
and the constant coefficients of the polynomials are then selected as the 
augmented states. An example 3rd order polynomial y=oco + ctix + a 2 x 2 + osx* 
is shown by the curved trajectory in Figure 7 for the approximation of a linear 
spring with backlash according to the piecewise linear trajectory. 

The polynomial approach for approximating the nonlinear function has several 

drawbacks as listed below. 

• The polynomial is a continuous approximation and it is difficult to 
estimate the size of the backlash from the coefficients ceo, cci, a 2 , a 3 . 
15 .The polynomial approach introduces unnecessarily additional augmented 
states and increases the complexity of the problem. 




One significant advantage of the augmented Kalman filter is the fact that the 
parameters to be estimated do not have to be constants. The polynomial 
20 coefficients are replaced by a single augmented state for the spring and 
backlash force x k . The augmented filter then estimates the force tk as a 
function of time and no particular shape of the underlying nonlinearity has to 
be assumed. 
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The complete augmented model can be written as follows. 

X\ = X3+V1 
X2 = X4+V2 

X 3 = — (U - D{X 3 - X 4 ) - p - Tf - Kr) + V3 

x\ = — {D{x 3 -x 4 ) +p) + u 4 

p = V 5 

where variables Vj are noise processes acting on the state and parameter 
variables. The covariance matrix for the augmented noise process v is given as 

E(v(t)v T (t)) = 

where Qxx is the covariance of the noise on the state derivatives, while Q PP is 
the covariance of the noise on the augmented state derivative, ie. the derivative 
of the backlash and spring force. Qpp can be seen as the main tuning 
parameter of the augmented filter. If Q PP is set to zero, the parameter p will 
remain a constant equal to its initial value. The larger the value of Qp P , the 
faster the augmented state will be updated. 

List of designations 

A Arm 

D damper coefficient 

-cf motor friction torque 

K r force from P-controller on the motor position 

M Motor 

u additional torque input 

tk nonlinear backlash and spring function 



0 



0 
Qpp 
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